Enhanced stochastic oscillations in autocatalytic reactions 
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We study a simplified scheme of k coupled autocatalytic reactions, previously introduced by To- 
gashi and Kaneko. The role of stochastic fluctuations is elucidated through the use of the van 
Kampen system-size expansion and the results compared with direct stochastic simulations. Regu- 
lar temporal oscillations are predicted to occur for the concentration of the various chemical con- 
stituents, with an enhanced amplitude resulting from a resonance which is induced by the intrinsic 
graininess of the system. The associated power spectra are determined and have a different form 
depending on the number of chemical constituents, k. We make detailed comparisons in the two 
cases = 4 and = 8. Agreement between the theoretical and numerical results for the power 
spectrum are good in both cases. The resulting spectrum is especially interesting in the k = 8 
system, since it has two peaks, which the system-size expansion is still able to reproduce accurately. 
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I. INTRODUCTION 

Autocatalytic reactions have long fascinated physicists 
and chemists because of their unique features [l| . A 
chemical reaction is called autocatalytic if one of the re- 
action products is itself a catalyst for the chemical reac- 
tion. Part of the reason for the interest in these types of 
reactions stems from the fact that even if only a small 
amount of the catalyst is present, the reaction may start 
off slowly, but will quickly speed up once more catalyst 
is produced. If the reactant is not replaced, the pro- 
cess will again slow down producing the typical sigmoid 
shape for the concentration of the product. All this is 
for a single chemical reaction, but of greater interest is 
the case of many chemical reactions, where one or more 
reactions produce a catalyst for some of the other reac- 
tions. Then the whole collection of constituents is called 
an autocatalytic set |2l ■ In addition to the interesting 
properties of autocatalytic sets, there is also an intrigu- 
ing possibility that "bootstrap" reactions such as this 
may have had an important role in producing complex 
or self-replicating molecules required for the origin of life 
on Earth 

Theoretical studies of the properties of autocatalytic 
reactions are typically of two kinds. In the first, rate 
equations for the reactions are written down and these 
are either solved numerically or their properties investi- 
gated using the techniques used in the study of dynam- 
ical systems. An alternative is to carry out computer 
simulations of the actual reactions themselves. However 



there is a third possibility: using methods from the the- 
ory of stochastic processes an analytic approach to the 
full model (and not just the mean field version) is possi- 
ble. In the last few years this last approach has been used 
for systems which are closely related to autocatalytic re- 
actions, such as predator-prey interactions 0], metabolic 
reactions [1], and epidemic models These all show 
oscillatory behavior in the number of individuals or con- 
stituents, which arise from feedbacks. These oscillations 
are distinct from the limit cycles found in the rate equa- 
tions, and are purely stochastic in origin. The main tool 
that is used to analyze these systems is the system-size 
expansion of van Kampen |Tl[ which gives very good 
agreement with the simulation results, even for systems 
of a moderate size. 

In this paper we apply this technique to the autocat- 
alytic reaction scheme studied by Togashi and Kaneko 
[Ta [Tsj . In most autocatalytic reactions there are two 
types of constituent: the autocatalytic and the substrate. 
The number of the latter type are kept constant by con- 
tinually feeding them in, however the former are not in- 
jected nor extracted from the system. In this sense the 
system is closed as far as the autocatalytic constituents 
are concerned, but open for the substrate. In the scheme 
that Togashi and Kaneko investigate, the reactions are 
cyclic, with k constituents Xi, . . . , Xk reacting accord- 
ing to X^ + Xi+i 2Xi+i with Xk+i = Xi, i = 1,. ..,k. 
The chemicals are assumed to be in a container which is 
well-stirred, but with the possibility of diffusing across 
the surface of the container into a particle reservoir. 
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In their approach Togashi and Kaneko [T3, IIjI use only 
computer simulation to study this reaction scheme. The 
analytic techniques we will use begin by writing down 
the master equation for their reaction scheme, and then 
studying it through a systematic expansion in A^~^/^, 
where N is the system size. To leading order one finds 
the rate equations which appear in (l2| . and to next-to- 
leading order a Langevin equation which describes the 
fluctuations about the stable fixed point of the rate equa- 
tions. From previous work we expect that (i) this first- 
order correction will be sufficient to yield results which 
are in good agreement with simulation data, (ii) the large 
amplitude of the oscillations can be understood as a 
resonant effect. One of the strengths of the technique 
is that the next-to-leading order corrections give linear 
Langevin equations which can be analyzed exactly for 
arbitrary values of k. 

The outline of the paper is as follows. In Section |TT] 
we derive the equations which govern the dynamics of 
the system, both in the deterministic limit and for the 
fluctuations about this limit. These fluctuations are an- 
alyzed in Section [TTTl bv calculating the power spectra for 
each chemical species i. Theoretical predictions are then 
compared to direct simulations for the case k — A and 
k = 8. Finally in Section IIVI we sum up and discuss 
possible future work. An Appendix contains the inter- 
mediate steps required to find the equations given in the 
main text. 



II. GOVERNING EQUATIONS 

The autocatalytic reaction scheme described in Section 
U can be formulated as 



E % X,; 



Ti + l 



2Xi+i , Xk+1 = Xi 



X,, 



E, i = l,. 



(1) 



Here rj,ai and /3i (with rk+i = ri), are the rates at 
which the reactions take place and E is the null con- 
stituent. Such constituents have to be included so that 
the number of molecules of type Xi, Ui, are all indepen- 
dent. If the size of the system is denoted by N , then 
Si=i ''^i + "--E — N , where is the number of null con- 
stituents. While N is fixed, ue may vary as the total 
number of molecules changes with time. In practice, ue 
does not explicitly appear in the formalism; it is always 
replaced by TV — 2_/i=i '^j- The rate constants and j3i 
in Eq. ([T]) represent the interactions of the system with 
the particle reservoir outside the container. In effect ai 
and (3i are the rate at which molecules appear and dis- 
appear from the system, and thus are analogous to birth 
and death rates. 

As an aside, we note that reaction rates which result 
from a binary encounter should be scaled by the vol- 
ume of the system, V . That is, Vi Vi/V . This follows 
from a straightforward kinetic theory argument [14i] . This 



is an innocent modification as far as this study is con- 
cerned, since it is carried out at constant volume, but 
it becomes crucially important when the volume is al- 
lowed to change, as it does in the analysis of the phase 
transition reported in [l^, [3 ■ 

The state of the system is labeled by the set of in- 
tegers {ni, . . . , rife} and, under the assumption that the 
transitions from this state to any other only depends on 
these integers, the system is Markov and may be de- 
scribed in terms of a master equation. In constructing 
the master equation we need to give the transition rates 
T{n'\n) from the state n to the to the state n', where 
n = (ni, . . . ,rtfc). If the system is well-stirred, so that 
the probability of a reaction taking place is proportional 
to its rate and the number of reactant molecules, then 
from Eq. ([1]) these transition rates are 



r(ni 



l,ni+i -I- 1, . . .,nk\n) = r^+i 



N N 



T(ni, . . . ,ni + \, . . . ,nk\n) = a 
T{ni, ...,ni-l,.. .,nk\n) = 



N 



(2) 



The master equation for the probability that the system 
is in state n at time t, Pin, t), may now be written down: 



dP{n,t) 
di 



- E - 1) 



X [T{ni, ...,ni- + 1, . . . ,nk\n)P{n,t)] 

k 

+ J2 (^r' - 1) [T{ni, . . . , n, + 1, . . . , nk\n)P{n, t)] 

i=l 
k 

+ ^ - 1) [T(ni, . . . , - 1, . . . , TZfe |n)P(n, t)] (3) 



where Sf^^ are the step-operators introduced by van 
Kampen p^ : 

ffV(n) = /(ni,...,n,±l,...,nfc). (4) 

Equations such as ([3]) are difficult to analyze, but if 
one is particularly interested in large or moderately sized 
values of N , then the system-size expansion provides an 
elegant way of encapsulating the essential aspects of the 
model. The key assumption of the method is to write 

m 



N 



N 



(5) 



From this relation, limN^cxiini /N) — (pi{t), the fraction 
of the molecules which are of type Xi at time within the 
mean- field {N oo) limit. The fluctuations about these 
are assumed to be Gaussian, hence the 1/\/N in Eq. 
One of the consequences of this assumption is that one is 
looking at a regime sufficiently far from boundaries that 



3 



the probability density functions of the Xi are Gaussian. 
This imphes that stochastic extinctions will not be well- 
described by the method, at least to leading order. 

Substituting Eq. ([5]) into Eq. ^ allows us to expand 
the master equation as a power series in 1/ ^/N. To see 
this we first note that the step operators ffl) take a par- 
ticularly simple form within the method [lO| 



1 d 



(6) 



If we set P{n, t) equal to n(^, t), the left-hand side of the 
master equation becomes ,10j 



dP{n,t) _ dU{tt) 



dt 



dt 



dnj^, t) del), 

d^^ dt 



(7) 



Substituting Eq. ^ into the right-hand side of the mas- 
ter equation ([3]), and using the transition rates ([2]), we 
may equate terms of the same order in 1/\/N on the 
left- and right-hand sides. To leading order this gives 
(see Appendix \K\ for details) 

d(j)- ( \ 

= [ri4>i_i - ri+i(/ii+i) -I- Qfi I 1 - ^ j - Pijii , 

(8) 

where r is a rescaled time: r — t/N . At next order one 
finds a Langevin equation: 



i=i 



(9) 



where M is a fc x fc matrix which can be found from 
Eqs. (|A5|) and (jA7p . and rji is a Gaussian white noise 
with zero mean and correlator 



(10) 



and Bij is another k x k matrix given by Eq. (|A6p . 

The first equation, Eq. ([81) , is a deterministic equation 
for the fraction of molecules which are of type i. It agrees 
with that of Togashi and Kaneko [l2|, if one takes into 
account that their equations are for concentrations and 
so contain the (constant) concentrations of the species in 
the reservoir. There is also an additional term (j)j in 
Eq. ([8]), which is typically present when mean- field equa- 
tions are derived in systems with a fixed size, but not in 
the phenomenologically postulated form. For small con- 
centrations it will not be important, but clearly it will 
have an effect as the ceiling on particle numbers is felt, 
reducing the number of molecules entering the container 
from the reservoir, as it should. The second equation, 
Eq. ([HI), is a stochastic differential equation for the devi- 
ation from these mean-field values. It is the analysis of 
these two equations that allow us to describe the stochas- 
tic aspects of the autocatalytic reactions in a quantitative 
way. 



III. ANALYSIS OF THE FLUCTUATIONS 

In their numerical studies, Togashi and Kaneko (l^.[l3j 
looked at the simplest case of the model where the rates 
TijOfi and j3i were the same for all chemical species. To 
illustrate our method we will do the same, and so from 
now on we will drop the subscript i on these constants, 
but it should be clear that our analysis also applies to the 
general situation where they are different for each species. 
With this choice, the deterministic equations ([8]) have a 
single fixed point: 



a 



f3 + kc 



(11) 



where the asterisk denotes the fixed point value. 

If N is so large that the fluctuations are completely 
negligible, then the system tends towards a state where 
the fractions of the chemical species in the system are 
equal, and given by Eq. pT|) . and stays there. Of course, 
if N is finite this is no longer the case and there are 
fluctuations about this stationary state — and as we will 
see these can be significant even if N is quite large. Since 
these fluctuations are expected to be oscillatory, we begin 
their analysis by taking the Fourier transform of Eq. ^ 
to find 



k 

E 



(12) 



where the / denotes the Fourier transform of the function 



/. Defining the matrix 
solution to Eq. is 



Mij to be ^ij{uj), the 



|,(c.) = ^$r.i(a;)^^.(w). 



(13) 



To identify the dominant frequency of the oscillatory 
behavior, we compute the power spectrum for the ith 
species, Pi{(-o), from Eq. p3|) : 



^.(^)^ 1^(^)1 



k k 



3 = 1 1=1 



(14) 

Since $ = —icul — M, where / is the k x k unit matrix, 
and since M and B are independent of the structure 
of Pi{uj) is that of a polynomial of order 2k divided by 
another polynomial of degree 2k. The explicit form of 
the denominator is | det $(aj)p. 

From previous investigations of fluctuations of a simi- 
lar kind 0, H, HI , we expect that the fluctuations about 
the stationary state pT|) will be enhanced by a resonant 
effect: for values of uj for which | det *&(w) | is a minimum, 
the power spectra will show peaks which correspond to 
larger than expected fluctuations at that frequency. This 
effect was first conjectured by Bartlett [l^ in the context 
of the modeling of measles epidemics, and later elabo- 
rated upon by Nisbet and Gurney [lB|, who called these 
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stochastically induced cycles, quasi-cycles. However it 
is only in the last few years that explicit calculations 
within the system-size expansion have been carried out 
and a quantitative understanding of the phenomenon has 
emerged 0- 

To understand the analytic structure of the power 
spectra, we begin by supposing that we can neglect the ef- 
fects of the numerator on the right-hand side of Eq. 
and simply determine the dominant frequency by look- 
ing for the value which minimizes | det $(w)|. The effect 
of the numerator will be to shift this frequency; we are 
assuming as a first approximation that this shift will be 
small, as indeed it has been found to be in some cases [71 . 
If Xj are the eigenvalues of M, then the denominator of 
the expression for the power spectra may be written as 

k 

\ det ^uj)\^ ^Y[{-iuj~- \j){iuj- X*) . (15) 

Since M is real, the Xj will be real or come in complex 
conjugate pairs, so that the products in Eq. psp has one 
of two forms: 

(i) If A is real, the two factors involving this eigenvalue 
give (w^ + A^). 

(ii) If A is complex: A ~ Xr + zA/, the four terms in- 
volving A and A* give 

\uj^ + {Xl-Xj)+2iXRXi\\ (16) 

The resonant effect has its origin in the structure of the 
factor coming from the complex eigenvalues shown in the 
expression ([1^ . It is smallest, and so gives the largest 
contribution when it is in the denominator, for frequen- 
cies which satisfy 

^c=A|-A|. (17) 

If there are several pairs of complex eigenvalues and their 
conjugates, the largest contribution should come from the 
pair for which XrXi is smallest. Clearly this will only 
be approximately true since, not only are we neglecting 
the numerator, but also the factors (w^ -I- A^) coming 
from real eigenvalues, as well as those coming from other 
complex conjugate pairs. However, as we will now see 
by looking at two specific cases, k = 4 and k = 8, these 
approximations appear to be remarkably good. 

We study the cases fc = 4 and fc = 8 because they 
are the smallest even values of k for which one complex 
conjugate pair and two distinct complex conjugate pairs, 
respectively, exist (there are two complex conjugate pairs 
for A; = 6, but they are equal, and three for k — 8, but 
two of these are equal). We therefore expect to see one 
peak in the power spectra when fc = 4 and two when 
fc = 8. Our analysis, and the accuracy of our approxi- 
mations, can be directly checked by numerical simulation 
of the chemical reaction system ([T]) by use of the Gille- 
spie algorithm [3, [l3| ■ This produces realizations of the 



stochastic dynamics which are equivalent to those found 
from the master equation ([3]). Averaging over many of 
these realizations gives us power spectra after Fourier 
transformation, which are exact to a given numerical ac- 
curacy. Wc now investigate the two cases fc = 4 and 
fc = 8 in more detail. 



A. Power spectra when k — 4 

The time evolution of the species is depicted Figure [T] 
This clearly displays large oscillations which we aim to 
investigate analytically. Before beginning this analysis, 
we observe that species 1, 3 (odd) and 2, 4 (even) are 
paired together and move up and down from the refer- 
ence mean-field level in a synchronized fashion. This fact 
was already recognized in [l^ . [isj and shown to drive suc- 
cessive switches between the 1-3 or 2-4 rich states, close 
to the absorbing boundary, i.e. when a small number of 
molecules is simulated. The rate at which the changes oc- 
cur is controlled by the diffusion parameter. However, the 
details of the transitions stem from a purely dynamical 
effect which cannot be captured within the perturbative 
analysis developed here. 

Let us now turn to analytically characterizing the 
aforementioned oscillatory regime. To this end we begin 
by determining the eigenvalues of the M matrix. From 
Eq. (|A14|) . these are 

Ao = /3 4a , A2 = /3 , 

Ai = (3 + 2ir(j)*, Ag-A^. (18) 

Within the approximations we have discussed, we would 
expect that there should be a single peak in the power 
spectrum for any one of the chemical species at a fre- 
quency given by (see Eq. pT|) ') 

u;l = 4r^rf-P' = -^^-P'. (19) 

{f3 + 4a) 

In Fig. [2] we show the power spectrum (for the chemical 
species i = 2) found by averaging over 500 realizations 
from the Gillespie algorithm, together with that found 
from Eq. using the matrices B and M given in the 
Appendix. The good agreement between the simulation 
results and those found from applying the system-size 
expansion, shows that the method works well for — 
5000. The parameters used in this case were r = 10 
and a = (3 = 1/64, which gives a value of ~ 4 from 
Eq. From Fig. [5] we see this is a surprising good 

estimate for the position of the peak, given the significant 
frequency dependence which we have neglected to obtain 
the estimate (|17p . 

Another check of the accuracy of these approximations, 
and so of Eq. ([T7]) , is to imagine increasing the parameter 
(3 at fixed r and a, and asking when uJ^ will become 
zero, and so at what frequency will the peak in the power 
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spectra disappear. From Eq. ([19)) we estimate this to be 



2ra 

/3 ~ — or /3 

P 



(20) 



which equals 0.56 for the values of r and a used in Fig. [51 
Once again this agrees well with the full spectrum which 
predicts the peak to disappear at about the same value. 
As a final check, we measure the position of the peak 
from a set of simulations run at different values of r. Di- 
rect measurements (symbols) are compared to the theory 
(solid line) in Figure[31and are in good quantitative agree- 
ment. Again, we recall that adjusting the rate r can be 
equivalently seen as modifying the volume of the system, 
which is the setting investigated in [H, • 



B. Power spectra when fc = 8 

From Eq. ([A14[) . the eigenvalues of the M matrix are 

Ao = /3 + 8a , \i= 13, 

Ai = /3 + V2ir0* , At = AJ , 

\2 = P + '2ir(t>* , As AJ , 

A3 = /? + V2ird)* , As = A^ (21) 

Since there are two distinct complex conjugate pairs we 
would expect to find two peaks in the power spectra, one 
at ojI = 2r'^{(j)*f - and the other at uj^ = 4r2(0*)2 - 
/3^. For small /3, one peak will be at a frequency -\/2 
times the other. We would also expect that the peak at 




80 ' 

Time 

FIG. 1: (Color online) Time evolution of selected species, 
i = 1,2,3,4 in clockwise order for the case k = 4. Here 
r = 10, a = /3 = 1/64, iV = 8192. The dashed line indicates 
the mean field solution. The species display a clear oscillatory 
trend about their mean field values. A paired synchronization, 
(1,3) vs. (2,4) rich states, is also visible, as already observed 

in [11, [a. 




FIG. 2: (Color online) Power spectrum of species i — 2 when 
fc = 4. The analytical curve is shown as a solid line and 
the simulation (average over 500 independent realizations) as 
symbols. Here r = 10, a = /? = 10/64, iV = 5000. 




FIG. 3: (Color online) The position of the peak in the power 
spectrum (for species i = 2 when k = 4) plotted as function 
of the rate constant r. Symbols refers to the stochastic sim- 
ulations, while the solid line shows the analytical prediction. 
Here a = /? = 1/64, = 5000 



lower frequency would be larger than the one at higher 
frequency, since XrXj is smaller for the former. That is, 
the pole in the power spectra in the complex frequency 
squared plane is nearer to the real axis for the peak at 
lower frequency, and so should have a bigger effect. So, 
in summary, our approximations indicate that the peaks 
in the power spectra should be given by 



f3' 



4r^a^ 



^c2 



(/? + 8a)' 



P\ (22) 
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FIG. 4: (Color online) Power spectrum of the time series for 
species i = 2 when fe = 8. The analytical result (solid line) is 
superimposed onto the simulations (symbols), averaged over 
500 independent realizations. Here r = 200, a = 1.9, /3 = 2, 
iV = 7000. 



with the peak at w = Wci larger than the one at w = 
LiJc2 ■ The results of plotting the full spectrum found from 
Eq. and simulation results are shown in Fig. 2] for 
r = 200, a = 1.9 and (3 — 2. This corresponds to peaks 
ai uo — 31.2 and lo = 44.14, according to Eqs. (P^ . which 
once again agrees very well the the results displayed in 
the figure, as does the prediction that the peak nearest 
the origin should be the largest. 



IV. CONCLUSION 

Auto-catalytic networks are central in many different 
contexts and play an important role in intracellular bio- 
chemical reaction schemes. In this latter scenario, species 
are confined in a closed volume, delimited by the cellular 
membrane. Low concentration can occasionally develop 
resulting from the complex mutual interaction between 
microscopic actors. Under such conditions, fluctuations 
matter and the effects of the intrinsic discreteness need 
to be properly accounted for. In other words, continu- 
ous kinetic equations prove inadequate, finite size cor- 
rections becoming significant. These aspects were nu- 
merically substantiated by Togashi and Kaneko [l^, 
within the framework of a simplified system of k coupled 
autocatalytic reactions. 

In this paper we have taken this forward by studying 
analytically the associated master equation via a system- 
atic expansion in power of N~^/'^ , where N is the system 
size. To leading order, the mean-field rate equations are 
recovered, while higher order corrections enable us to ex- 
plain the large amplitude of the oscillations as detected in 
direct simulations. Importantly, the calculation applies 
to arbitrary values of k. For k — A a. peak in the power 
spectrum is found, while for fc = 8 two peaks develop. 
To the best of our knowledge, this is the first time that 
a double-peaked power spectrum has been predicted to 
emerge as a resonant effect, within a van Kampen type 



of analysis. In both cases, theory and simulations agree 
well thus confirming the importance of finite N contri- 
butions. Possible extensions of the present work include 
taking spatial variations into account. This could yield 
spatial oscillations in the species concentration, which 
would again be driven by the discreteness of the system 
components. 
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APPENDIX A: THE FINITE iV EXPANSION 

In this appendix we will give the intermediate steps of 
the calculation using the system size expansion, starting 
with the master equation ^ and ending with the results 
(I5))- ((T(I| . We will also give the explicit expressions for the 
matrices M and B. 

Applying the ansatz ([5]) to the right-hand side of 
Eq. ([3]), the step-operators ([U take the form the 
rii in the transition rates ^ are replaced by (pi and S^i 
using Eq. ^ and P{n,t) becomes n(^,t). This yields 
the following terms: 

(a) Terms of order N^-^^^: 

k 



i=l 



d 



d 



(b) Terms of order N ^ and involving first order 
derivatives: 

( d d d 

{ ri+i<l>i-^ii+i + Ti+icpi+i—ii ~ ri+i(l)i— — ^i+i 

d ^ d 9 1 

-ri+icj)i+i^ + tti V -rrp^] + Pi-^^i ) n(^, t) ■ 

Oc,i+i ~[ c'^i I 



(A2) 



(c) Terms of order N ^ and involving second order 
derivatives: 



- 2- 



92 



V 02 



-ai 



i-E*^^ ] ^ + ^^'^«^^n(|,i). (A3) 
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Introducing r = t/N^ the terms of order iV~^/^ in 
Eq. (|A1[) may be identified with the second term on the 
right-hand side of Eq. ^ . This gives the N deterministic 
equations ([5]). The terms of order in Eqs. (|A2[) and 
(|A3|) . are now identified with the remaining term on the 
right-hand side of Eq. This resulting equation is a 
Fokker-Planck equation: 



an 

97 



(A4) 



From Eq. (jA2p we see that the Ai(^) are hnear functions 
of the and from Eq. (jA3[) that the are independent 
of them. Exphcitly: 

^i(^) = {nfpi-i - ri+i(f>i+i) + ri(l)i^i^i 

k 

- r^+l(p^^^+l - ^ - , (A5) 



and 



Bii — 



a j = i-l 



ri+i4>i4>i+i + ri4>i4>i_i 



if j = i + 1 



(A6) 

In Eqs. (|A5p and (|A6p . (/>fc+i = (/>i and ^fc+i = ^i, which 
follows from the cyclic nature of the model. 

Since the Ai{^) are linear fimctions of the we may 
write them as 



MO 



(A7) 



This means that the probability distribution at next-to- 
leading order, n(^,T), is completely determined by the 
two k X k matrices Af and B, whose elements are in- 
dependent of the , and only functions of the (jfj . For 
our purposes, where we need to Fourier analyze the fluc- 
tuations, it is more convenient not to use the formula- 
tion in which the fluctuations are described by a Fokker- 
Planck equation, but rather in terms of Langevin equa- 
tions. The Fokker-Planck equation (|A4p is completely 
equivalent to the Langevin equation ^ with the corre- 
lator dioi) [Ulii, and it is this formalism that we will 
use. 

In principle the matrices M and B are time dependent, 
since 0^ is. However, in practice we are interested in fluc- 
tuations about the stationary state, and so we are only 
interested in the values that the (pj take on at late times. 



Furthermore, in Section IIIII we studied the simple case 
ri — r,ai — a and (3i — f3, for which the relevant value 
of the is given by Eq. (fTTj) . With these assumptions 
M and B are given by 



(A8) 
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7713 


77io 


777,1 


7772 • ■ 


. 7712 


7712 




7712 


7713 


7770 


7711 • ■ 


. 7712 


7712 


M = 


7712 


77i2 


7773 


7710 • ■ 


. 7772 


7712 




7712 


7712 


7712 


77l2 . . 


• mo 


7711 




7711 


7712 


7772 


7772 ■ ■ 


■ TO3 


771o_ 



where 

Too = - 
and 



-a— (3 , TOi — —a—rcj)* , 7712 = —a , TO3 = 



-a+r(j) , 
(A9) 



B 



bo h 0. 
bi bo 61 . 
bi bo. 



bi 





0. 
bi 0. 



bo bi 
bi bo. 



(AlO) 



where 



60 = 2r{(b*r + + a(l - kcj)*) , &i = )' • (All) 



We note that M is a circulant matrix [20|, and there- 
fore its eigenvalues are given by 



k 



A£ = > ;77iy e(""'(^-l)^)/^ ^ = 0,1,...,A:-1, (A12) 



where mij is the element of M in the first row and jth 
column. In fact, M is not the most general form of cir- 
culant matrix; (fc — 3) entries in each row are equal (to 
7712). This leads to a simplified form for the eigenvalues: 



A, = Too + 7711 e^"^/^- + TO3 + 7772 ^^^'''^^ 

i=2 

= TOO + 7.1 e--/^ + TO3 e---/^ - TO2 5^^^^^ , 



(A13) 



where in the last line £ ^ 0. Putting in the values from 
Eq. (M\i gives 



P + ka, if £ = 

(3 + 2ir(j)*sin{2Tre/k), if ^ 7^ . 



(A14) 
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